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Summary 

A lifting-surface theory has been developed for a 
helicopter rotor in forward flight for compressible and 
incompressible flow. The method utilizes the concept 
of the linearized acceleration potential and makes 
use of the doublet lattice procedure. Calculations 
demonstrating the application of the method are 
given in terms of the lift distribution on a one-bladed 
rotor, a two-bladed rotor, and a rotor with swept- 
forward and swept-back tips. Also, the lift on a rotor 
vibrating in a pitching mode at 4 per revolution is 
given. Compressibility effects and interference effects 
for a two-bladed rotor are discussed. 

Introduction 

The aerodynamics of rotating elements, such as 
helicopter rotors and propellers, has been under ex- 
tensive study since the advent of the airplane, and, 
with a combination of experimental and analytical 
approaches, succcessful designs have been achieved. 
In many cases, two-dimensional theory has been 
used, usually modified by assumed spanwise distri- 
bution and inflow velocities. This report presents a 
compressible, unsteady lifting-surface method for a 
helicopter rotor in forward flight within the limits of 
linear aerodynamic theory. 

The method is based on the concept of the ac- 
celeration potential, originally formulated for three- 
dimensional oscillating lifting surfaces by Kussner 
(1941). The method was first applied to an oscil- 
lating wing in uniform rectilinear motion, including 
effects of compressible flow by Runyan and Wool- 
ston (1957). The acceleration-potential approach has 
now become standard in the determination of the un- 
steady aerodynamic forces for flutter studies of lifting 
surfaces in rectilinear motion. 

An early use of the acceleration-potential ap- 
proach for a rotating system was made in a paper 
by Hanaoka (1962) for the loading on a marine pro- 
peller in incompressible flow. The acceleration po- 
tential has been used in the past for studying the 
propeller noise problem. In all these noise propaga- 
tion cases, the problem was specialized early in the 
analytical development to the so-called far-field case 
usually with a stationary observer. However, the 
lifting-surface theory is essentially concerned with 
the details of the near-field case for a co-moving ob- 
server and with the satisfaction of certain boundary 
conditions. Runyan (1973) utilized the acceleration- 
potential approach to obtain a solution for the os- 
cillating propeller in compressible flow. Suciu et al. 
(1976) have derived an incompressible lifting-surface 
theory and applied it to a windmill. The procedure 
is based on the velocity-potential method and subdi- 


vides the integration areas into panels within which 
a constant pressure distribution is assumed. Dat 
(1973) derived a general expression for an acceler- 
ation doublet for any motion. The theory developed 
by Dat was applied to a helicopter in forward flight 
by Costes (1973). The approach of Costes is similar 
to the method given in this paper, with the excep- 
tion of a numerical differentiation procedure which 
is adopted to obtain the downwash velocity from the 
velocity potential. In the present approach, the ex- 
pression for downwash velocity is obtained analyti- 
cally and includes a number of implicit differentia- 
tions. Pierce and Vaidyanathan (1983) treated the 
helicopter rotor in forward flight using the method 
of matched asymptotic expansion for the incompress- 
ible case. 

The method presented in this report sets forth 
a formulation of fundamental, three-dimensional, 
compressible, unsteady aerodynamic theory for pro- 
pellers and helicopter rotors in forward flight. An 
inertial coordinate system is adopted, and the in- 
tegrations involved in solving the integral equation 
are formulated for arbitrary space-time variations. 
In this formulation, singular terms arise in the in- 
tegrations and are handled by the technique of the 
finite-part integrals. 

The paper is divided into three basic sections. 
The first section contains a brief derivation of the 
fundamental equations and is followed by a section 
describing a method of solving the integral equation. 
Finally, the method is applied to specific examples, 
and the computational results are given. 

Symbols 

A, B coefficient in equations (28) 

A! rotor-blade area 

A nm aerodynamic-influence coefficients 

A n ,B n Fourier coefficients 

C chord of rotor 

Ct thrust coefficient (Thrust/7rpQ 2 i?^) 

where thrust is normal to tip-path 
plane and is positive in positive 
2 -direction 

c speed of sound 

D absolute value of D 

D radius vector from doublet to down- 

wash point 

D = D/D , unit vector of D 



h flapping deflection, positive up 

h 0 magnitude of flapping deflection 

I value of singular integral 

K kernel function 

L.E. leading edge 

M t i p tip Mach number 

£,m,n direction cosines of n 

l 0 ,m 0 , n a direction cosines of n 0 

n unit vector at downwash point, normal 

to velocity vector 

n 0 unit vector at doublet point, normal to 

velocity vector surface 

p pressure 

q source or doublet strength 

R s rotor root radius 

Rt rotor tip radius 

r distance of downwash point along span 

rg, r u lower and upper bounds of panel in 
spanwise direction 

r 0 distance of doublet along span 

r 0 distance of doublet along span at 

singular point time 

r 0 - (rg + r u )/2 

T.E. trailing edge 

t field time 

U velocity of rotor system, parallel to 

z-axis, positive in negative i-direction 

— * 

V velocity at downwash points 

V n resultant velocity normal to blade 

V 0 velocity of doublet 

W velocity of rotor system, parallel to 

2 -axis 

w n downwash velocity, normal to velocity 

plane, positive down 

X position vector of downwash point 

from inertial frame origin 

X 0 position vector of doublet from inertial 

frame origin 

x, y , 2 Cartesian coordinates of downwash 

point 


x 0 ,y 0 , z 0 Cartesian coordinates of doublet 
position 

x a distance from pitch axis to downwash 

point 

a twist angle at downwash point as a 

function of t 

o-o twist angle at doublet position as a 

function of r 

a r angle of axis of rotation relative to 

2 -axis 

l V 0 /c 

9 angular position of blade at time 
0 = Q t 

&B blade pitch angle relative to plane of 

rotation 

9 0 angular position of blade at time r, 

9 0 = Or 

9 n angle of velocity vector V n 

9 W blade angle of attack 

p advance ratio, U/QRt 

1 0 air density 

r,r 0 source time 

t source time at which integrand in 

equation (25) becomes singular 

d velocity potential 

velocity potential of a moving doublet 

<f> s velocity potential for a moving source 

'f r acceleration potential 

acceleration potential for a moving 
source 

d) azimuth angle measured from down- 

wind position 

fi rotation speed of rotor 

cj/ ( flapping frequency 

u> a pitching frequency 

A dot over a symbol indicates derivative with 
respect to time. 

Basic Formulation 

The formulation of the aerodynamic equations 
is based on the linearized acceleration-potential ap- 
proach. The fluid is considered perfect, with no 
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separation, and the formulation is based upon the 
assumption of small perturbations. One reason for 
adopting the acceleration-potential approach is that 
the pressure discontinuity occurs only on the sur- 
face of the blade, and thus the boundary conditions 
need only be applied on the blade surface, and not 
throughout the wake. The blade is treated as a very 
thin surface of discontinuity across which a pressure 
jump occurs. The effect of compressibility is taken 
into account by utilizing the complete linearized po- 
tential for a lifting doublet and by utilizing the effects 
of retarded time. 

As shown in sketch 1, an inertial coordinate sys- 
tem was used. The helicopter rotor moves in the neg- 
ative x-direction with velocity U, moves in the pos- 
itive ^-direction with velocity W, and rotates coun- 
terclockwise with a constant angular velocity SI. A 
point of interest on the rotor blade is designated by 
the vector X 0 (t) from the origin of the ground-based 
coordinate system. 


use the relationship between the velocity potential 
and the acceleration potential to provide an equation 
for applying the boundary condition of “no flow” 
through the rotor surface. The relationship between 
the pressure and velocity potential for an inertial 
coordinate system is 



( 2 ) 


where ^ is the substantive derivative. Dropping out 
the second-order terms and integrating with respect 
to observer time results in 

m = f* ( 3 ) 

J — OO 


Taking the derivative of <j> with respect to a di- 
rection n gives the fluid velocity component in the 
direction ft as follows: 


Z 



Sketch 1 

If $ is the acceleration potential of a source (or 
doublet), the perturbation pressure is then given by 

p=-p9 (1) 

This expression represents the pressure p at point 
X due to a single source (or doublet) located at 
X 0 . The potential contains a constant value of q, 
which represents the strength of the source and thus 
the magnitude of the pressure. In this form, there 
is no boundary condition available to determine the 
value of the arbitrary constant q and the resulting 
pressure. Since the spatial derivative of a velocity 
potential represents a velocity, it is convenient to 


A w n ~^r = jr[ * dt ' ( 4a ) 

The right-hand side of equation (4a) represents an 
induced velocity in the direction n due to a moving 
source (or doublet) of strength q where the value of 
q is contained in the expression for By adding 
the contribution of all doublets distributed over the 
rotor, the induced velocity along n at the downwash 
point can be represented as 

Wji — j Aw n dA' = j ~-f* $ dt' dA' (4b) 

The left-hand side, w n , represents the known bound- 
ary condition and is the velocity normal to the veloc- 
ity vector at the downwash point. (A discussion of 
w n is included in the section entitled “Definition of 
Boundary Conditions.”) Thus, the problem resolves 
itself by setting up a method of solution of equa- 
tion (4b) from which the values of q, the unknown, 
can be determined which satisfy the known velocity 
boundary conditions w n . The next section contains 
the development of the expression for an acceleration 
doublet potential. 

Acceleration Potential of a Doublet in 
Compressible Flow 

The acceleration potential $ s must satisfy the 
wave equation 

= /(*,*) (5) 
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where f{X,t) is a source distribution. Furthermore, 
if the path of an isolated source is a variable of a 
function of time X 0 (t), then f(X,t) — q6[X — X 0 (t)}, 
where 6 is the Dirac delta function. 

Using the Green’s function formulation, the 
acceleration-potential expression for a moving source, 
'Fs can be written as (Morse and Feshbach 1978, 
page 841) 


Equation (9) was also derived by Dat (1973) in 
a different fashion. For incompressible flow, c — ► 
oo, the first term — > 0, and the integral remains 
unchanged except for the upper limit, where r = t. 

To obtain the final equation for downwash Aiu n , 
a second directional derivative is required. This 
derivative is taken normal to the flight path at the 
location of the downwash point as follows: 




Q(Xo,T) 


An \X - X 0 (t ) j 


! V 0 (tHX-X 0 (t)1 

c]x-X 0 (t) 


(6) 


Aw n 


<Hp 

dn 


= ft • Vx<t>D 


This leads to 


where X 0 {t ) designates the position of the source at 
time r, X is the position of the field point at time 
t, V 0 (t) is the velocity of the source point at time 
r, c is the speed of sound, and q is the strength of 
the source. An auxiliary equation which relates the 
time interval (t — r) to the distance between the two 
points is given as follows: 


t — t = 




( 7 ) 


Aw n — 


(H • n Q — n Q ■ Dn * D) 


4ircD 2 (X-b>j$) 
n 0 - fin • D - ® • h 0 n • D - n 0 ■ Dn • D -f- H 0 • Dn • fi 

(i - 6 • J) 

n 0 ■ bn ■ 6(i - p 2 + & ■ 0) 


+ 


(1- D-0) 2 
r T o(r 0 ) 


Ho • Dn * Dq 


T 0 


47tc 2 P(1 -D‘fi ) 2 


To 


± f To[ro) (n-fio bfi‘ b\ 

4 *J-oo q \ ° 3 ) 


dr 


( 10 ) 


and is usually referred to as the causality condition. 
Equation (6) expresses the potential as a function 
of r, and only through equation (7) as an implicit 
function of t and X. 

From equations (3) and (6), the velocity potential 
due to a moving source is 


f at'= r 

J — OO J —{ 

=-/ T 

An 

J —oc 




4n(D — D • 0) 


dt' 


Qifi) 

D(r') 


(8) 


where D = X-X 0 , D = |D| , and dt' = (l - dr'. 

The quantities t 1 , t', and t, r satisfy equation (7). 

A source potential cannot be used to produce a 
pressure difference across a lifting surface. However, 
a doublet does contain the proper form of the singu- 
larity on the surface to provide a jump in pressure 
over the lifting surface. The expression for a doublet 
potential can be obtained by taking the derivative of 
a source potential in a direction normal to the air- 
foil surface. If n 0 designates the normal to the airfoil 
surface, then 

Q 

<t>D M = 4>s {t) = n 0 * = —Ho • sy^4>s 

= — ^ 4 - 1 + r q ^idr' ( 9 ) 

Anc(D - p ■ D)D\ t y_ 00 4 ttD3 


This concludes the derivation of the basic equa- 
tion for the downwash Aw n at a field point created 
by an arbitrarily moving doublet. The term contain- 
ing q in the present application was not used in sub- 
sequent calculations, since the rate of change of the 
loading is relatively small. Consequently, this term 
compared with the other terms in the brackets is rel- 
atively small. It can be argued that the ratio of the 
q term to the fifth term in the first bracket is QD/c, 
which is small unless the frequency Q is very large, 
which may be the case in other applications, such 
as acoustics. The remaining portion of this paper is 
directed towards specializing this equation for a heli- 
copter blade moving with a forward velocity U, and 
a vertical velocity W, and rotating with an angular 
velocity fi. 

Equation (10) gives the downwash at a field 
point ( x,y,z,t ) due to a doublet placed at a point 
(x 0 . Vo, z 0 , t) with a strength q. To form a lifting sur- 
face such as a rotor, it is necessary to distribute the 
doublets over the lifting surface and integrate over 
the surface to obtain the downwash at a field point. 
If the downwash is known, the quantity q can then 
be determined. Letting K be the expression on the 
right-hand side of equation (10), the final equation is 



( 11 ) 
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where A' is the area of the rotor surface. This equa- 
tion represents a rather formidable computing task, 
and the history of lifting-surface theory, even for non- 
rotating wings, has centered on devising approximate 
methods to accomplish the integration in an econom- 
ical manner. One method, termed the doublet lattice 
method, has been very successful when applied to air- 
craft wings and is probably the most economical pro- 
cedure of the many variants. This method was first 
demonstrated for the unsteady case by Runyan and 
Woolston (1957) and was later expanded by Albano 
and Rodden (1969). This is the method adopted in 
this paper, and the application is discussed subse- 
quently. In the next section, the coordinate system 
is specified, and the doublets and downwash points 
are appropriately located for a helicopter rotor in for- 
ward flight. 

Specification of Coordinate System 

The blade has the chord C and length Rt~ R s , 
where R s is the radial distance to the root of the 
blade and Rt is the distance to the tip of the blade. 
Sketch 2 shows the blade momentarily coinciding 
with the coordinate system along the positive z-axis 
at t = 0 and at a new location at time t later. 
The blade executes a counterclockwise rotation 



with angular velocity fi while moving with velocity 
U along the negative ^-direction and velocity W 
along the positive 2 -direction. Since the doublet 
lattice method has been adopted, the doublet point 
for one chordwise panel lies 1/4 C ahead of the 
midchord, and the downwash point lies 1/4 C aft 
of the midchord. The positions of the doublet point 
and downwash point can be established as follows. 
The Cartesian components of the doublet position 


are 


x 0 — —Ur + r 0 cos(Qr) — {C/4) sm(Qr) cos a 0 
yo — To sin(nr) + (C/4) cos(Or) cos a 0 > 

z 0 == Wt + {C/4) sina 0 


( 12 ) 


where r 0 is the radial distance of the doublet along 
the span. (See sketch 1.) With the substitution of 
C — > —C,r 0 — ► r, and r — > t, the position of the 
downwash point is given by 


x = —U t + r cos(fit) + (C/4) sin(fit) cos a ' 
y = r sin(fit) — (C/4) cos(fit) cos a > 

z = Wt — (C/4) sin a 


(13) 


In equations (12) and (13), the angles a and a 0 are 
the twist angles of the velocity vectors V and V 0 , 
respectively, defined by 


W 'j 

ana U sin(f2f) -I- rO . . 

w r (14) 

tan (Xq r 

C sin(f2r) -I- r 0 fl J 

The reference plane defined by the doublets and 
downwash points is a twisted surface. The doublet 
velocity, namely the time derivative of the position 
vectors, can be computed from equation (12) as 
follows: 

V 0 = x 0 t + y 0 j + z 0 k 

The unit vector n 0 is chosen to be perpendicular 
to the twisted surface (see sketch 3) created by the 
velocity V 0 , which is a function of r 0 (see eqs. (12) 
and (14)). To determine n 0 as a function of r, write 



Sketch 3 

Vo = V' 0 + Wk, where V' 0 is the projection of the 
doublet velocity on the horizontal plane. Then, 

n 0 — t 0 i + m 0 j + n 0 k (15) 
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where £ 0) w 0 , and n 0 are the directional cosines of 
the unit vector n 0 . Form the vector quantity S by 
making £? = x W k. Therefore, 


§ = W[U 4 ■ r 0 0 sin(nr) ■+ (<7/4)0 cos(Ot) cos a 0 ]1 

4 W[rotl cos(Or) — ((70/4) sin(Or) cos a]i (16) 


Furthermore, m 0 , and n G can be obtained by 
imposing the following conditions: 


n 0 LS ' 
ft 0 ±V 0 > 

Wo\ = 1 i 


(17) 


From equations (15), (16), and (17), 


£o — 


W[U 4 roO sin (fir) 4 (<7/4)0 cos (Or) cos a 0 ] ) 

VoV^ + Vo' 2 


W^[r 0 0 eos(Or) — (<7/4)0 sin(Or) cos a 0 ] 
m 0 — — 


Vq\/w^ +VP 


n 0 — 


K 


yfw* + VP 

where V 0 is the magnitude of V 0 , and 

V/ 2 — [17 + r 0 nsm(nr) + (C , /4)f2cos(nr)cosa 0 ] 2 

+ [r 0 ncos(fJr) - (C/4)fisin(nr)cos£* 0 ] 2 (19) 


By the same procedure, ft is derived by making the 
appropriate substitutions, C — *• — C, a 0 ot,r — >■ t, 

etc. (i.e., n = li + rrij + nfc), to get 


f _ W[U + rfisin(nt) — (C/4)0cos(nt) cos a] ' 

Vi\Zw*+V* 

_ — WjrO cos(Ot) + (<7/4) 0 sin (Ot) cos a] 

v 0 v^ 2 +Vo 4 

>/iy2 + y'2 


( 20 ) 


where 


y /2 — [</ 4 rO sin(Ot) — ((7/4)0 cos(Ot) cos a] 2 
4 [rO cos(Ot) 4 (<7 /4) 0 sin(Ot) cos a ] 2 

The vector D = X — X 0 defined in equation (8) can 
be expressed as 

D = ({LT(t - r) 4 rcos(Ot) - r 0 cos(Or) 

4 (C/4) [sin(Ot) cos a 4 sin(Or) cos a 0 ]} 2 
4 (rsin(Oi) — sin(Or) 


■— (C/4) [cos(Oi) cos a 4 cos(Or) cosa 0 ]} 2 
4 [W(t - r) - (C/4) (sin a 0 4 sina)] 2 ) 1 ^ 2 (21) 

With the substitution of the quantities n 0 , n, F 0 , V, D, 
and D into equation (11), the integral equation was 
solved for the unknown <?(r 0 , r) by using a collocation 
process based on the doublet lattice assumption. The 
kernel is singular when D == 0, and this was handled 
by use of the finite-part technique. 

Solution of Integral Equation 

The basic integral equation to be solved 
(eq. (11)) is 

Wn — j I K dA' 

A' 


This equation states that doublets of strength q are 
distributed over the area of the rotor surface A 1 , and 
the induced velocity normal to the blade surface is 
equated to the normal velocity of the blade, where 
w n is the boundary condition on the blade surface 
and is a known quantity. Thus, the only unknown in 
the integral equation is the strength of the doublet q. 
As mentioned previously, one method for obtaining 
a solution to this equation is termed the doublet 
lattice method and is used in this report. It has been 
used successfully for the wing in pure translational 
motion. One of the advantages of the doublet lattice 
method is that it eliminates the chordwise integration 
and thus reduces computing costs. The reduction is 
the result of a judicious choice of the location of the 
doublet and the downwash points. 

Doublet Lattice Technique 

In following the doublet lattice technique, the ro- 
tor is divided into a number of panels, both span- 
wise and chordwise. In each panel, a line of doublets 
of unknown strength q^ is located at the 25-percent 
chordwise location of the particular panel, and the 
downwash is evaluated at the 75-percent chordwise 
location and midspan of the panel (see sketch 4) for 
a rotor divided into 2 chordwise panels. The problem 
is then resolved by performing an integration from 
to r u for the vortex at C/8 for the downwash point in 
the middle of the spanwise panel at 3C/8 (w n = uq). 
The effect of the same vortex must be determined for 
the downwash located at 7C/8 (u> n — u^)- A simi- 
lar calculation is made for all the remaining vortices, 
including lateral distances. Therefore, a collocation 
procedure is used to obtain a set of equations in terms 
of the unknown loadings q{. 
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Sketch 4 


The term q(r 0 , r) represents the strength of the 
doublet located at r 0 and at time r, and is propor- 
tional to the unknown loading. Within a particular 
panel, q is assumed to be constant in the spanwise di- 
rection. To account for unsteadiness, a solution was 
formulated to take into account the time variation of 
the strength of the doublet. The time variation of 
q is represented by assuming a Fourier series of the 
form 


m 

g(r 0 ,r) — A 0 (r 0 ) + [An (r 0 ) cos(nn-r) + B n (r 0 ) sin(nfir)] 
n=l 

( 22 ) 

If q(r 0 ,r ) is assumed to be a function of r 0 alone, 
which means that the doublet strength does not vary 
with time, the Fourier series reduces to q{r 0 ) = 
A 0 . A solution obtained with this approximation is 
termed the quasi-steady solution. 

This series (eq. (22)) was inserted into the basic 
equation (eq. (11)) and integrated with respect to r. 
However, there were more unknowns than simulta- 
neous equations to solve for the unknowns. The ad- 
ditional required equations were obtained by evaluat- 
ing equation (11) at a number of azimuth locations. 
For instance, if m — 1, then 

q(r 0 , t) = A 0 + A\ cos(fir) 4- B\ sin(fiT) (23) 

The azimuth was divided into equal segments of 120°, 
and the proper boundary conditions were applied 
at ip = 0°,120°, and 240°. Thus, the necessary 
additional equations were provided. 

A set of equations is obtained as shown below: 


Wl = Anqi + Ai242 H 1- MmQm \ 

w 2 — ^2191 + A 2292 H + M mQm 


(24) 


W n — A n iqi + A n 2Q2 + * * ' + A nm q m j 


where A nm = / r r “ K nm dr Q and where n refers to the 
downwash point and m refers to the panel number. 

The rotor can be divided into any number of 
segments both chordwise and spanwise, with more 


segments resulting in more accurate results. How- 
ever, the computing time and costs increase rapidly 
(roughly proportional to the square of the numbers 
of segments) as the number of segments is increased; 
therefore, a balance must be considered in any given 
situation, and there must be a trade-off between com- 
puting costs and desired accuracy. 

Numerical Integration of Kernel 

A closed-form integration of the kernel has not 
been found; therefore, the integration was performed 
by numerical integration, except for the area sur- 
rounding the singularity. The integration domain 
was divided into five areas. (See sketch 5.) The in- 
tegral over areas 1 to 4 (hatched) were computed 
numerically by using the two-dimensional Romberg 
quadrature (Davis and Rabinowitz, 1984) and the 
contribution of the singular region (unhatched) was 
obtained in closed form by consideration of the finite 
part, as shown in the next section. 



The integral in the downwash equation (eq. (10)), 
is singular when D — * 0 and produces a complication 
which must be properly treated. The integration 
along r is the path the doublet has taken in arriving 
at the final doublet point at (C /4, r 0 ) measured in the 
local blade coordinates. The integration takes place 
along the path from r = -00 to the final doublet 
position at r 0 . The distance D is the distance from 
the integration point at time r to the downwash point 
at X, 
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Sketch 6 


There is a particular set of values of r 0 and r for 
which the denominator D approaches zero, and the 
result is an infinite integrand. The function is given 
above in sketch 6. The singular part of the kernel is 

I = jf "" jf* * ' n ° - 3 <^ iTo dr „ (25) 


A :== [U (t — t) + r cos# + C/4(sin 0 cos a + sin #<> cos <* 0 )] 2 
4 - [rsin# — C/4(cos#coso: + cos# 0 cos a 0 )} 2 
-f [W(t — t) — C/4(sina 4- sina 0 )] 2 (28b) 


In this equation, a 0i defined in equation (14), is a 
function of r 0 . By using the doublet lattice method, 
the rotor is divided into spanwise segments from rg 
to r u . If these segments are small, a 0 will vary 
approximately as 


da 0 ^ WQ 

dro ** ~(Usm0 o + r o n ) 2 + W 2 


(29) 


If the value of a 0 is approximated as a constant and 
evaluated at midsegment f 0 , it is possible to integrate 
equation (26) in closed form in the r 0 direction. This 
is quite acceptable in the helicopter mode. That is, 
both a 0 and da 0 /dr 0 are small. Therefore, reversing 
the order of integration gives 



n • n 0 

(r 2 — 2 Br 0 + A) 3 / 2 


dr 0 dr 


(30) 


As D — » 0 at the downwash point, D becomes per- 
pendicular to n; therefore, at the singular point, the 
second term is zero and is neglected in the treat- 
ment of the singularity. However, this second term 
is retained in all the numerical integrations involving 
areas 1 to 4 because it represents an important inter- 
ference term, particularly when the blade is passing 
over a trailing wake. The time and distance at which 
I becomes singular are designated by t and r 0 . The 
domain of the integration in equation (25) consists 
of a rectangle in which the duration r 2 — t \ is kept 
extremely small. In other words, the integration is 
performed along a slit in r 0 , over which the second 
term in equation (25) is negligible. Therefore, the 
integral I can be approximated by 

'-nr**** <*> 

By squaring the quantities in D (eq. (21)) and 
collecting coefficients of r 0 , the distance D can be 
written as 


D = (r%-2Br 0 + A) 1 / 2 (27) 


The vector n 0 is also a function of r 0 , but in the re- 
gion of the singularity it has a very small variation 
and is evaluated at the singular position (f,f 0 ). Per- 
forming the r 0 integration results in 

j _ f T 2 n • n 0 ( ru + fi/2 
Jti A — B 2 ^ y/rQ — 2Br u + A 

UMH 1 *„ (si) 

yj r t 2 4 - 2 Brg 4- A J 


The singular term in this expression is A - B 2 . 
Equation (31) may also be written as 


= P 7 n* 

Jn f{r) 


where /(r) is the singular part and /(r) = A — B 2 
and g{r) is the nonsingular part as follows: 


/(t) = (?7(t - r) — r sin [fi(t — t)\ + C/ 4 {cos [f 1(t — t)] cos a 
2 

4- cos olq} ) 4- [W (t — t) — C/Asm a + sin a 0 ] 2 (32) 


where 

B = cos 0 o [U{t — r) 4- r cos# 

+ C/4(sin#cosa + sin#o cosao)] + sin#o[rsin# 

- C/4(cos#cosa 4- cos# 0 cosa^)] (28a) 


Equation (32) shows that /(t) is nonnegative, so the 
zero of /(r), at r = f, is a second-order zero. 

Expanding /(r) in a Taylor series about the sin- 
gular point f yields 

/(r) = /(f)+/'(f)(r-f)+/"(f)(r-f) 2 /2+- • ■ (33) 
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Since f(f) is a second-order zero, 


Definition of Boundary Conditions 


/(f) = /'(f) =0 (34) 

The conditions given by equation (34) have been 
verified numerically. If only the square term is kept 
in equation (33), equation (31) can be written as 


1 = 



2 

1W) 


qjr) g'(f) g"(f) 

(r - f) 2 (r - f ) 2 J 


(35) 

In equation (35), if t<i and ri are chosen symmetri- 
cally about f , the odd derivative terms integrate to 
zero. Futhermore, the third term can be neglected, 
since g n {f) is small. The major contribution comes 
from the first term. Then, using the standard inte- 
gration technique (Mangier 1951), the final result for 
the integral is 


r _ g(0 4 

/"(f) Ar 


(36) 


where 2 At = T 2 — ri, and t\ < f < T 2 - 

A numerical problem arises, because this finite- 
part integration results in a negative number close 
to the total of the surrounding numerical integration 
areas, which are positive. Thus, it is necessary 
to take the difference between large numbers, and 
the final integration accuracy is dependent on the 
accuracy of the two integrations. On the one hand, 
the numerical integration is more accurate if Ar 
is kept large, because the very large values of the 
integrand near the singularity are avoided. On the 
other hand, regarding the finite-part integration, the 
denominator was expanded in a Taylor series about 
the singular point f. Therefore, it is desirable to 
maintain At as small as possible to keep within the 
limits of the applicability of the series expansion. 
Numerous calculations were made, varying At until 
a reasonable convergence was found. The value used 
for the calculations in this report was 0.01(t — f), 
or 1 percent of the time difference. Actually, there 
is little difference between 1 percent or 10 percent of 
the time difference, and the computing time and cost 
are considerably reduced by using 10 percent. For 
trend studies 10 percent is recommended principally 
to reduce computer costs. However, for final-design- 
type analysis, 1 percent is recommended. 

For the spanwise direction, A r 0 is also an 
integration-limit variable. The finite-part integral 
was obtained by approximating the angle of twist of 
the velocity vector across a segment. This approxi- 
mation is made by assuming that the angle of twist 
is constant across the segment, with a value deter- 
mined at the center of segment. Numerical experi- 
mentation indicates that for a helicopter, Ar 0 = 0 is 
satisfactory. 


The integral equation (eq. (11)) equates the in- 
duced velocity caused by a distribution of doublets 
(right-hand side) to the velocity normal to the ve- 
locity vector (left-hand side) to form the boundary 
condition of no flow through the blade. The preced- 
ing sections have shown the development of the aero- 
dynamic induced velocity in terms of the unknown 
loadings q^. This section contains a definition of w n 
for several conditions. 



Sketch 7 

Sketch 7 illustrates the geometric condition at a 
particular azimuth ip and radial distance r on the 
blade. The angle 0g represents the physical angle 
of the blade relative to the x-y plane. The angle 6 n 
represents the angle of the velocity vector V n normal 
to the centerline of the blade. The angle 0 W is the 
blade angle of attack; therefore, 

Ow = ~~ @n ( 3 7 ) 

where 

_i W 

0 n = tan . . (38) 

(rU + U sm Ot) v ’ 

The velocity normal to the velocity vector V n is 

Wn,o = Vn tan Oyj (39) 


where 

V n = y/^n + U sin 0 f) 2 + W* 

If the blade has a built-in twist, then 8q is a function 
of r and would be calculated at the appropriate 
downwash points. 

For a rotor having a swash plate, 


0 S — Aq-Aj cos(nt)-Bi sin(0£)— cos(2nt)-B2 sin(20£)H 

(40) 
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An additional contribution to the normal velocity 
results from the rate of pitching of the blade. If the 
distance from the pitch axis to the downwash point 
is denoted by x a , the normal velocity to be added is 

W n ,s = x a 0 s (41) 

For a wing harmonically vibrating in a pitching 
mode defined by 0 a , 

Ob = 0 a cos(w a t + <f>a) (42) 

where 0 a is the deflection mode shape and is a func- 
tion of r, oj a is the frequency of pitching oscillations, 
and <fi a is the phase angle relative to ^ = 0. The 
additional velocity to be added to w n is. 

w n,a = %oi'0b (43) 

For a rotor oscillating harmonically in a flapping 
mode given by 

h = h 0 cos (v/jf 4- (fin) (44) 

where h 0 is the magnitude of the spanwise flapping 
deflection, which is a function of r, where ui^ is the 
flapping frequency, and where <fih is the phase angle 
relative to \l> .= 0, 

w n,h = h (45) 

Therefore, for a rotor having a swash plate de- 
fined by equation (40), vibrating in a pitching mode 
(eq. (42)), and vibrating in a flapping and pitching 
mode (eq. (44)), the final downwash equation would 
be 

Wn = w n ,o + w n>8 + w n>a + w n ^ h (46) 

Also, the entry of the rotor into a gust field could be 
handled or a propeller rotating close to the fuselage 
could be treated if the velocity field through which 
the rotating surface were known. 

Application to Specific Examples 

The foregoing analysis has been applied to sev- 
eral specific examples. (See fig. 1.) The following 
section presents results for several paneling config- 
urations (e.g., 5 spanwise panels and 1 chordwise 
panel, designated (5-1) and 7 spanwise panels and 3 
chordwise panels, designated (7-3)). The untwisted 
rotor blade was maintained at a constant pitch set- 
ting $b of 0.1 rad for all the calculations. Although 
the spacing of the paneling can be arbitrary, the pan- 
eling arrangements used in all examples in this paper 


are uniform except for the swept-tip cases, for which 
the arrangements are described. 

Incompressible Flow 

Single blade. To investigate the convergence of 
the method of this report when using the doublet 
lattice procedure, the program was run for several 
chordwise and spanwise elements for the incompress- 
ible case. The thrust coefficient Ct is plotted against 
the azimuth angle in figure 2. (In all the plots for 
thrust coefficient versus azimuth angle, the thrust 
was calculated for 16 uniformly spaced azimuth an- 
gles, and each curve was faired using a cubic spline.) 
The rotor was first divided into 5 spanwise panels 
and one chordwise panel (5-1), and the results are 
shown by the solid line. The chordwise division was 
increased to (5-2), and the results are shown by the 
long dashed line. It can be seen that very little 
change occurred. The spanwise divisions were in- 
creased to (7-1), and the largest change between var- 
ious panelings occurred at xfi = 0°, where the dif- 
ference in Cj is about 11 percent. Increasing the 
chordwise divisions to 3 (7-3) produced little change 
from the (7-1) case. 

An interesting phenomenon occurs in the rotor 
first quadrant. For ip = 0 to 37°, the lift increases to 
a local maximum at ip = 37°, abruptly falls to a local 
minimum for %fi = 60°, and then rapidly increases to 
a maximum at tfi = 100°. A similar phenomenon is 
shown analytically by Egolf and Landgrebe (1983) in 
figure 60 of that report, where a local minimum and 
a local maximum occur in the same range of azimuth 
angles, even though the geometry of the two blades 
and the flight conditions are different. In figure 93 
of the same report, some test data show a similar 
variation of loading in the same azimuth range. 

The chordwise pressure distributions for the (7-3) 
case are presented in figure 3. In using the doublet 
lattice method, the loading is concentrated at the 
location of the doublet which for the (7-3) case is 
located at 0.0833(7, 0.416(7, and 0.75(7. The pressure 
was faired using a cubic spline through the three 
vortex locations and the known value of zero at 
the trailing edge. The distributions are given for 
7 spanwise positions. In general, the curves exhibit 
the expected shape, with the largest values occurring 
as the leading edge is approached. For the span 
distribution, the values at r/Rt — 0.85 are slightly 
larger than the values at r/Rt = 0.95. These values 
indicate a falling off in the tip region. 

From these concentrated forces, the section pitch- 
ing moment can be calculated. Figure 4 presents 
these results for ip = 90°. The section moment was 
taken about the 1/4 C and a nose-down moment is 
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taken as positive. The pitching moment shows some 
rather dramatic changes along the span. The mo- 
ment is nose-up near the tip ( r/Rt = 0.95), changes 
to a small nose-down value, then becomes nose-up for 
most of the inboard region. Integration of the mo- 
ment would result in a total root-end pitching mo- 
ment up at ip = 90°. 

Swept tip . The segments used for the doublet 
lattice for the swept tip studies were (5-1), where two 
equal segments were used in the tip region and three 
equal segments were used in the unswept inboard 
section. Lift is plotted against azimuth in figure 5 
for the two sweep conditions and for zero sweep. In 
general, the three results show little difference. The 
swept-back configuration has the largest lift from 
ip k 300° to 40°. For ip * 100° to 240°, the 
swept-forward configuration has a slightly larger lift 
than the other configurations. The total lift for 
one rotation for the swept-back case and the swept- 
forward case appear to give about the same lift as 
produced by the unswept rotor. In figure 6, the lift 
distribution along the rotor span is given for ip = 0°. 
The major effect of sweep is concentrated at the tip, 
where the swept-back tip load is greater than either 
the unswept or swept-forward cases. In contrast, 
figure 7 shows a larger forward tip load than for the 
unswept or swept-back tip for ip = 180°. 

Blade oscillating in pitch . An example of un- 
steady loads on a rotor blade with (5-1) paneling 
oscillating in a pitching mode about the midchord 
at a frequency of 4 per revolution (120 cycles/sec) is 
given in figure 8. For this case, a 17- term Fourier 
series (m = 8) was used to simulate the oscillating 
load, which was comprised of one constant term, 8 co- 
sine terms, and 8 sine terms. The nonoscillating and 
oscillating rotor-blade loading is given for one revo- 
lution. The blade was oscillated through an angle of 
±0.1 rad about a mean angle of 0.1 rad. The effect 
of the oscillation is readily apparent when compared 
with the nonoscillatory case. With the harmbnic rep- 
resentation of the loading, the magnitude and phase 
of the several harmonic loads are easily determined. 
The magnitudes are plotted in figure 9. The only 
harmonic loads that were significantly changed from 
the steady case were the 3rd, 4th, and 5th. All three 
harmonics were increased, but the 4th harmonic was 
drastically increased. Another calculation was made 
for the nonoscillating case and was compared with 
the quasi-steady case. Virtually no difference was 
observed; this similarity indicates that, at least for 
this case, the rate of change of loading in a revolu- 
tion of the blade is small enough that the effect of a 
variable wake is negligible. 


Compressible Flow 

One-blade rotor . For a one-bladed rotor (5-1), 
the effect of compressibility is illustrated in figure 10, 
in which C? is plotted against azimuth angle. The 
incompressible result is included for comparison. As 
expected, the compressible load is larger than the 
incompressible load throughout one revolution and 
effect is greatest in the region of the advancing blade 
and smallest in the retreating region. 

Two-bladed rotor. The method has been extended 
to the two-bladed rotor (5-1 per blade) for the com- 
pressible case, and the results are shown in figure 11. 
The thrust coefficient C t per blade is plotted against 
azimuth angle for a one-bladed rotor and for a two- 
bladed rotor. For azimuth angles from ip = 20° to 
120°, the one-bladed rotor has a larger C For 
ip == 120° to 260°, Ct on the one- and two-bladed 
rotors are approximately the same. However, for 
ip = 260° to 340°, a drastic reduction in lift per blade 
occurs for the two-bladed rotor compared with the re- 
sults for the one-bladed rotor. The lowest lift occurs 
at ip = 292°; this places the other blade of the two- 
bladed rotor at ip = 112°, the point of maximum lift. 
Apparently, the high lift on the blade at ip = 112° 
creates a very unfavorable induced velocity on the 
second blade at ip = 292°, which requires the load- 
ing to go to zero to satisfy the boundary conditions 
at ip — 292°. 

Concluding Remarks 

A linearized lifting-surface theory, including the 
effects of compressibility, has been developed for a 
helicopter rotor in forward flight. The method uti- 
lizes the concept of the acceleration potential and 
makes use of the doublet lattice procedure for per- 
forming the required integrations. Also, the method 
has been extended to include the effects of unsteady 
flow and blade motion. 

Sample calculations have been done for several 
cases. These include the effect of swept-back and 
swept-forward tips. The effect of these two tip config- 
urations was minimal on the total loading for one rev- 
olution. However, the loading distribution changed 
considerably for several azimuth positions ip . Com- 
pressibility was investigated for one configuration. 
As expected, the effect was greatest in the advancing- 
blade region (ip = 90°) and was minimal in the 
retreating-blade region. A comparison of the thrust 
coefficient Cp of a one-bladed rotor and a two-bladed 
rotor was made. In the azimuthal range between 
20° and 120°, the one-bladed rotor showed higher 
lift. However, between ip = 260° and 340°, the two- 
bladed rotor showed lower Cy. The effect on C? of 
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a blade oscillating in pitch at 4 per revolution has 
been shown. The effect on the total blade lift has 
been shown and the effect of the oscillation is readily 
apparent. The harmonic content was calculated, and 
the greatest difference between the oscillatory and 
nonoscillatory cases was in the 4th harmonic. 

NASA Langley Research Center 
Hampton, VA 23665-5225 
September 3, 1985 
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Figure 1. Configurations used in examples and input parameters. 
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Figure 2. Four panel configurations for a single rotor blade. Incompressible; p = 0.17; 0g — 0.1 rad; 
a r = 0.05 rad; O = 30 rad/sec. 
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Figure 3. Chordwise pressure distribution for several spanwise locations, ip = 90°; incompressible; p = 0.17; 
9b = 0.1 rad; a r = 0.05 rad; f) = 30 rad/sec. 



Figure 4. Spanwise section moment distribution about 1/4 C. Positive nose down; ip = 90°; incompressible; 
p = 0.17; 9b = 0.1 rad; a r = 0.05 rad; fi = 30 rad/sec. 
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Figure 5. Comparison of lift on a swept-back zero-sweep and swept-forward blade. Incompressible; n = 0.17; 
6b = 0.1 rad; a r = 0.05 rad; Q = 30 rad/sec. 



Span, r/R t 

Figure 6. Spanwise section lift distribution for swept-tip configurations, ip = 0°; incompressible; // = 0.17; 
Ob =0.1 rad; a r = 0.05 rad; O = 30 rad/sec. 
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Figure 7. Spanwise section lift distribution for swept-tip configurations, rp = 180°; incompressible; n — 0.17; 
Oq = 0.1 rad; a r = 0.05 rad; n = 30 rad/sec. 


X 10 3 



Azimuth angle, deg 

Figure 8. Comparison of lift on a rotor blade oscillating in pitch, at 4 per revolution with lift on a nonoscillating 
blade. Incompressible; fx = 0.17; = 0.1 rad; a r = 0.05 rad; 17 = 30 rad/sec. 
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Harmonic number 


Figure 9. Harmonic content for nonoscillating and oscillating cases. r/Rt= 0.95; incompressible; fi = 0.17; 
6q = 0.1 rad; a r = 0.05 rad; O = 30 rad/sec. 
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Figure 10. Incompressible and compressible lift for one-bladed rotor, /z = 0.17; 6b= 0.1 rad; a r = 0.05 rad; 
M tip = 0.54. 
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Figure 11. Lift on two-bladed and one-bladed rotor plotted against azimuth angle. Compressible; fx 
9g = 0.1 rad; a r = 0.05 rad; M t j p = 0.54. 


.8 



1. Report No. 2. Government Accession No. 

NASA TP-2503 

4. Title and Subtitle 

Compressible, Unsteady Lifting-Surface Theory for a Helicopter 
Rotor in Forward Flight 

7. Author(s) 

Harry L. Runyan and Hsiang Tai 


9. Performing Organization Name and Address 

NASA Langley Research Center 
Hampton, VA 23665-5225 


3. Recipient’s Catalog No. 

5. Report Date 

December 1985 

6. Performing Organization Code 

505-33-43-09 

8. Performing Organization Report No. 

L- 15976 
To. Work Unit No. 

11. Contract or Grant No. 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 


13. Type of Report and Period Covered 

Technical Paper 

14, Sponsoring Agency Code 


15. Supplementary Notes 


16, Abstract 

A lifting-surface theory has been developed for a helicopter rotor in forward flight for compressible and 
incompressible flow. The method utilizes the concept of the linearized acceleration potential and makes 
use of the doublet lattice procedure. Calculations demonstrating the application of the method are given 
in terms of the lift distribution on a one-bladed rotor, a two-bladed rotor, and a rotor with swept-forward 
and swept-back tips. Also, the lift on a rotor vibrating in a pitching mode at 4 per revolution is given. 
Compressibility effects and interference effects for a two-bladed rotor are discussed. 


17. Key Words (Suggested by Authors(s)) 

Helicopter 

Rotor 

Lifting surface 
Unsteady flow 
Compressible flow 
Acceleration potential 
19. Security Classif. (of this report) 

Unclassified 


18. Distribution Statement 

Unclassified — Unlimited 


Subject Category Ofx" 


20. Security Classif. (of this page) 

Unclassified 


21. No. of Pages 22. Price 

19 A02 


For sale by the National Technical Information Service, Springfield, Virginia 22161 


N AS A-Langley, 1985 



